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Abstract 

The relation between wildfire orientation and size is analyzed by means of a nonparametric 
test for directional-linear independence. The test statistic is designed for assessing the indepen- 
dence between two random variables of different nature, specifically directional (fire orientation, 
circular or spherical, as particular cases) and linear (fire size, scalar), based on a directional- 
linear nonparametric kernel density estimator. A bootstrap resampling procedure is provided, 
in order to apply the proposed methodology in practice. The performance of the test is checked 
through a simulation study and applied to analyze wildfire data from Portugal. 

Keywords: Bootstrap; Directional-linear density; Independence test; Nonparametric estimation; Wildfires 
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1 Introduction 

In the past decades, wildfires have been the major force of landscape change in Portugal, as in 
many other mediterranean type climates (Moreira et al. 2001; Silva et al. 2011). In Portugal, 
high primary productivity, dry and hot summers, abundant ignitions, and profound socio-economic 
changes characterized by rural abandonment and consequent homogenization of landscape with fire 
prone-shrublands, have combined to favor an increase in occurrence of large and devastating wild- 
fires (Moreira et al. 2011). Roughly 80% of annual burnt area is associated with very large fires 
occurring in 10% of days in the fire season, under the influence of well characterized atmospheric 
circulation pattern dominated by a ridge over the Iberian Peninsula (Pereira et al. 2005). Never- 
theless, Costa et al. (2011) showed that in Portugal, and over the period between 1980 and 2000, 
the relative importance of human and climatic drivers was temporal and spatially variable. Their 
results point towards the existence of a dynamic interaction between socioeconomic and landscape 
fire drivers, suggesting that while climate promotes favorable fire conditions across the entire coun- 
try, socioeconomic and landscape factors determine a large portion of the complex fire patterns 
found at regional scale. 

Independently from their drivers, wildfire patterns at landscape scale have important management 
implications (Moreira et al. 2001; Lloret et al. 2002; Moreira et al. 2011). For instance, it has been 
widely shown that landscape fuel reduction treatments will only be successful if strategically placed 
in order to intersect fire spread in the heading direction (Finney 2001; Schmidt et al. 2008). 

Barros et al. (2012) assessed the existence of preferential fire perimeter orientation at watershed level 
to support the spatial layout of fuelbreak networks. This analysis identified clusters of watersheds 
where fire perimeters are preferentially aligned along the NE/SW and the SE/NW axes. Those 
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watersheds include fire perimeters that together account for roughly 65% of overall burnt area in 
Portugal, over the period between 1975 and 2005, while in the remaining watersheds fire perimeters 
are aligned randomly. The authors argued that spatial patterns of fire perimeter orientation found 
in the 31-year dataset could be explained by dominant weather during the Portuguese fire season 
(Pereira et al. 2005). However, given that fire perimeter orientation analysis is event -based, i.e. it 
is based on the orientation of each fire event, all perimeters are treated equally independently of 
their size. In this paper, a test for assessing independence between wildfire size and orientation is 
presented, extending the work of Barros et al. (2012). Furthermore, orientation of the wildfire will be 
considered in a two-dimensional and a three-dimensional space. In Figure 1, some descriptive maps 
on the data of interest are displayed. The left plot shows the total area burnt in each watershed, 
whereas the middle plot represents the mean slope of the fires in each region. Finally, the right plot 
shows which watersheds exhibit a preferred fire orientation, versus a random orientation, according 
to Barros et al. (2012). 
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Figure 1: Descriptive maps of wildfires in Portugal with the 102 watersheds delimited by Barros et al. (2012). 
The left map reproduces the number of hectares burnt from fire perimeters associated with each watershed. 
Each fire perimeter is associated with the watershed that contains its centroid. The center map represents 
the mean slope of the fires of each watershed, where the slope is measured in degrees (0° aims for plain slope 
and 90° for a vertical one) . Finally, the right map shows watersheds where fires display preferential alignment 
according to Barros et al. (2012). 

Spatial characterization of a wildfire, by means of its main orientation, and the associated burnt 
area must be handled by non-standard statistical approaches, given the special nature of the fire 
orientation. Particularly, the wildfire orientation may be measured as an angle in the plane (two di- 
mensional orientation) or as a pair of angles identifying a direction in the three dimensional space. 
Hence, appropriate methods for handling circular and, more generally, directional data must be 
considered, jointly with suitable combinations of directional and linear techniques. 

The analysis of the relation between circular and linear variables has been classically approached 
through the construction of circular-linear correlation coefficients. The initial proposal by Mardia 
(1976) was further studied by Liddell and Ord (1978), who obtained its distribution under different 
hypothesis. From a different perspective, circular and linear variables can be also jointly modeled 
by the construction of circular-linear distributions. In this setting, Johnson and Wehrly (1978) 
introduced a method for deriving circular-linear densities with specified marginals. A new family 
of circular-linear distributions based on nonnegative trigonometric sums, adapting the method by 
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Johnson and Wehrly (1978), was proposed by Fernandez-Duran (2007), which proved to be more 
flexible. Recently, Garcfa-Portugues et al. (2012a) exploited the copula representation of the John- 
son and Wehrly (1978) family, allowing for a totally nonparametric estimator, which is applied to 
analyze SO2 concentrations and wind direction in a site close to a power plant. Nevertheless, the 
aforementioned methods are designed for the circular-linear case, whereas in our context, a more 
general tool is needed, provided that wildfire orientation may be reported in two or three dimensions. 

In this work, the assessment of the relation between a directional (circular, as a particular case) 
and a linear variable is approached through the construction of a formal test to check directional- 
linear independence. Inspired in the ideas by Rosenblatt (1975) for the linear setting (see also 
Rosenblatt and Wahlen (1992) and Ahmad and Li (1997)), the proposed test statistic is based on a 
nonparametric directional-linear kernel density estimator and an £2 distance is taken as a discrep- 
ancy measure between the joint estimator and the one obtained under the independence hypothesis. 

This paper is organized as follows. In Section 2, some background on kernel density estimation, for 
linear, directional and directional-linear data is presented. Section 3 is devoted to the introduction 
of the test statistic, describing in detail its practical application by means of a resampling procedure 
and a simplified version of the statistic, enabling its practical use. The finite sample performance 
of the test, in terms of size and power, is assessed through a simulation study on circular-linear 
and spherical-linear settings. Application to real data is provided in Section 4, including data 
description and results, focusing on the assessment of independence between wildfire orientation 
and burnt area size in Portugal. Some discussion and final comments are given in Section 5. 

2 Background 

In the linear setting, the basic building block for the independence test introduced by Rosenblatt 
(1975) is a kernel density estimator. Independence between two linear random variables is as- 
sessed through an £2 distance between a bidimensional kernel density estimator and the product 
of the kernel density estimators of the marginals. For the extension of such a procedure to the 
directional-linear case, kernel density estimation for linear, directional and directional-linear vari- 
ables is required. In this section, a brief background on kernel density estimation is provided. First, 
the linear kernel density estimation is introduced, providing some classical references on this topic. 
The extension of kernel density estimation for directional data is presented afterwards, following 
the proposals by Hall et al. (1987) and Bai et al. (1988). Finally, the directional-linear kernel 
density estimator, studied by Garcfa-Portugues et al. (2012b) is described. Some brief remarks on 
bandwidth selection are also included. 

2.1 Linear kernel density estimation 

The well-known kernel density estimator for linear data was introduced by Rosenblatt (1956) and 
Parzen (1962). From a random sample Zi, . . . , Z n drawn from a linear random variable Z (i.e. with 
support Supp(Z) C M) with density /, the kernel density estimator at a point z G R is defined as 



where K is a kernel function, usually a symmetric density about the origin, and g > is the 
smoothing or bandwidth parameter, which controls the roughness of the estimator. Properties 
of this estimator have been deeply studied (see Silverman (1986) or Wand and Jones (1995) for 
comprehensive reviews). It is also well known that, from a practical point of view, the considera- 
tion of different kernels (Gaussian, Epanechnikov, etc.) does not present an impact in the overall 
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shape of the estimator. However, the bandwidth is a key tuning parameter. Large values of the 
bandwidth parameter will produce oversmoothed estimates of /, whereas small values will provide 
undersmoothed curves. Thus, an appropriate choice of the bandwidth parameter is always an issue. 
In the testing problem that will be presented, this selection cannot be simply made by seeking the 
minimization of a certain error criterion for the estimator. Since several kernel estimators will be 
involved, along with a resampling procedure, further considerations on bandwidth selection will be 
detailed in the next section. 

2.2 Directional kernel density estimation 

The previous kernel density estimator has been adapted to different contexts, including directional 
data. Introduced by Hall et al. (1987) and Bai et al. (1988), following two different perspectives 
in the treatment of directional data, the directional kernel density estimator presents the same 
constructing elements as the classical one, that is, a kernel function and a bandwidth or smoothing 
parameter. First, denote by X a directional random variable with density /. The support of such 
a variable is the g-dimensional sphere, namely Q, q = {x G M. q+1 : x\ + ■ ■ ■ + Xq+\ = l}> endowed 
with the Lebesgue measure in Q q , that will be denoted by oj q . Therefore, a directional density is a 
nonnegative function that satisfies f n /(x) Wq(<ix) = 1. 

Given a random sample Xi, . . . , X n , of a directional variable X with Supp(X.) C Q q and density /, 
the directional kernel density estimator at a point x G 9 is given by 

i=l 

where L is the directional kernel, h > is the bandwidth parameter and Ch )q {L) is a normalizing 
constant depending on the kernel L, the bandwidth h and the sphere dimension q. The scalar 
product of two vectors, x and y, is denoted by x'y, where ' is the transpose operator. Note that the 
inclusion of the scalar product implies that distances between points in the sphere £l q are measured 
through the cosine of the angle between them. 

A particular feature about the kernel function is that directional kernels are not directional densities 
but functions of rapid decay and, in order to guarantee that the resulting estimator is indeed a 
directional density, the normalizing constant Ch, q (L) is required. In this sense, a common choice for 
the directional kernel is L(r) = e~ r , r > 0, also known as the von Mises kernel due to its relation 
with the q-von Mises-Fisher distribution (Watson 1983). In a (/-dimensional sphere, the von Mises 
density vM((i, n) is given by 

/„Af(x;/i, k) = C q (n) exp{Kx'/Lt}, 

(2) 

being fx £ £l q the directional mean, k > the concentration parameter around the mean and Z„ 
the modified Bessel function of order u, 

Then, if the kernel considered is von Mises (the negative exponential), the value of Ch, q (L) is 
C q (l//i 2 ) 1 e 1 /^ 2 and the directional estimator (1) can be interpreted as a mixture of q-von Mises- 
Fisher densities as follows: 

1 n 

A(x) = -^2 fvM (x; Xi, l/h 2 ) , 
i=i 
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where, for each von Mises component, the mean value is the i-th observation Xj and the common 
concentration parameter is given by l/h 2 , involving the smoothing parameter. Note that large 
values of h provide a small concentration parameter, which results in a uniform model in the sphere, 
whereas small values of h give high concentrations around the sample observations, providing an 
undersmoothed estimator. Cross-validation rules for bandwidth selection in the directional context 
were discussed by Hall et al. (1987), and for the particular case of circular data, plug-in rules were 
introduced by Taylor (2008) and Oliveira et al. (2012). 

2.3 Directional linear kernel density estimation. 

Consider a directional-linear random variable, (X, Z) with support Supp(X., Z) C £l q x M and joint 
density /. For the simple case of circular data (q = 1), the support of the variable is the cylinder 
and, in general, the support is a multidimensional cylinder. Following the ideas in the previous 
sections for the linear and directional cases, given a random sample (Xi, Z±) , . . . , (X n , Z n ), the 
directional-linear kernel density estimator at a point (x, z) G fl q x K can be defined as: 

t < \ c Kq (L) ^ Tv ( l-^ z-Zj \ 

A,,(x, z) = — g LK — ) , (3) 

where LK is a directional-linear kernel, g > is the bandwidth parameter for the linear component, 
h > is the bandwidth for the directional component and Ch, q (L) is the normalizing constant for the 
directional part. The previous estimator was introduced by Garcia- Portugues et al. (2012b), who 
also studied its asymptotic properties in terms of bias and variance, and established its asymptotic 
normality. 

A product kernel LK(-, •) = L(-) x if (■), specifically, the von Mises-normal kernel 

LK(r,t) = e' r x <p{t), r G [0, oo), t G R, 

will be considered along this paper in order to simplify computations, where (j> denotes the standard 
normal density, Af(0, 1). Similarly to the previous linear and directional kernel density estimators, 
a smoothing parameter (bidimensional, in this case) is involved in the estimator construction. The 
cross-validation procedures introduced by Hall et al. (1987) in the directional case can be adapted to 
the directional-linear setting, providing likelihood cross-validation better results than least squares 
methods (see Taylor (2008) or Oliveira et al. (2012) for some discussion for circular data). 

For an illustration of the kernel density estimators, a data sample from a circular-linear variable has 
been drawn. The linear component follows a mixture of normal distributions, whose kernel density 
estimator can be seen in Figure 2 (right panel). The circular component is given by a mixture of 
two von Mises and the circular kernel density estimator, with von Mises kernel, is plotted in Figure 

2 (central panel). Finally, contour plots for the circular-linear kernel density estimator can be seen 
in Figure 2 (left panel). Bandwidth parameters have been chosen by likelihood cross-validation 
criteria for the circular and the circular-linear case. A plug-in rule has been used for the linear 
kernel density estimator (Sheather and Jones 1991). Note that the support of the circular-linear 
variable is the cylinder, so the contour plot displayed in Figure 2 could be wrapped-around along 
the horizontal axis, with a perfect match of the left and right sides of the contour map. 

3 A test for directional— linear independence 

The test statistic for assessing independence between a directional and a linear variable will be 
described in this section. The proposal extends the idea by Rosenblatt (1975) to the directional- 
linear setting. Although the construction of the estimators is fairly simple, there are some issues 
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Figure 2: From left to right: kernel estimations for the circular-linear density and the corresponding circular 
and linear marginal densities, respectively. The circular-linear density is the three-mixture \vM ((1, 0), 2) x 
W(0, 0.5) + §uM ((0, 1), 10) xW(l, l) + \vM ((-1, 0), 2) xA r (2, 1), from where the circular and linear marginals 
are obtained. Solid lines aim for the true densities, where dashed lines for the kernel estimators. For the 
circular-linear and circular estimators, bandwidth is chosen by likelihood cross-validation, and for the linear 
estimator, by the plug-in rule of Sheather and Jones (1991). For the three densities, random samples of size 
n = 200 are drawn. 

that may hinder its application in real data analysis. First, the asymptotic distribution of the 
test statistic may not be useful from a practical perspective, given the slow rate of convergence, 
an issue which is quite common when dealing with tests based on smoothing methods. Hence, an 
alternative route relies on a proper design of a resampling scheme. For the problem considered, it 
will be seen that two difficulties arise: (i) the resampling mechanism, based on a smooth bootstrap 
approach, must regard for the balance between the resampling and the estimation bandwidths; (ii) 
an easy to compute expression for the test statistic must be available, overcoming the integration 
over multidimensional cylinders. Solutions to these two issues will be also provided, together with 
an illustration of the test performance by a simulation study. 

3.1 The test statistic 

Consider the joint directional-linear density f(x,z) f° r the variable (X, Z). Denote by /x and fz 
the directional and linear marginal densities, respectively. The null hypothesis of independence 
between both components can be stated as 

H : / ( x,z)(x,z) = /x(x)/z(z), V(x, z) Gl],xl 

and the alternative hypothesis as 

H a ■ /(x,z)(x, z) / /x(x)/z(^), for any (x, z) G Q, q x R. 

Following the idea of Rosenblatt (1975), a natural statistic to test Hq arises from considering the £2 
distance between the nonparametric estimation of the joint density f(x,z) by the directional-linear 
kernel estimator (3), denoted by f(x,z);h,gi an d the nonparametric estimation of f(x,z) under Hq, 
given by the product of the marginal directional and linear kernel estimators, denoted by fx-,h an d 
fz,g, respectively. Therefore, the following test statistic can be considered: 

T n = A 2 (f(X,Z);h,g, fx.;hfZ;g^) , (4) 

where A 2 stands for the squared £2 distance in £l q x R between two functions /1 and ji'- 

A 2 (/i,/ 2 )= / (/i(x,z)-/ 2 (x,z)) 2 ^(dx)dz. 
Jn q xR 
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The test statistic depends on a pair of bandwidths (h,g), which is used for the directional- 
linear estimator, and whose components are also considered for the marginal directional and lin- 
ear kernel density estimators. Under the null hypothesis of independence, Hq, it is clear that 

E [/(X 1 Z);ft lfl (x,z)] = E [/X;h(x)] E [f Z -g{z) . 

Asymptotic properties of (4) have been studied by Garcia- Portugues et al. (2013), who proved its 
asymptotic normality under Hq, but with a slow rate of convergence that does not encourage its use 
in practice. In addition, the construction of T n requires the calculation of an integral over £l q x R, 
which may pose computational problems related to the fact of solving several nested integrals. 
However, if the kernel estimators involved in the test statistic are obtained using von Mises and 
normal kernels, then an easy to compute expression for T n can be obtained, as stated in the following 
lemma. 

Lemma 1. If the kernel estimators in (4), obtained from a random sample {(Xj,Zj)}" =1 o/(X, Z), 
are constructed with von Mises and normal kernels, the following expression for T n holds: 

T n = l n (\*(h) o Sl(g) - \*(h)n(g) + \*(h)l' n l n n(g)) l' n , (5) 

where o denotes the Hadamard product and ^(h) and £l(g) are n x n matrices given by 

*W=[ n ^TZ\/h^ ^{g) = {^ g {Z l -Z ] )) ip 



being l n = (1, . . . , 1)', with length n, C q the normalizing function (2) and 4>^2 g the normal density 
with zero mean and standard deviation V2g. 

The proof of this result can be seen in Appendix A. Note that the previous formula for T n does 
not require integral computations, and it is only based in matrix operations. Hence, this will be 
the expression used for computing the test statistic within the resampling method for calibrating 
its distribution under Hq, as it provides a significative computational acceleration. 



3.2 Testing in practice 

Bootstrap calibration is not at all foreign to hypothesis testing, specially if the test statistics are 
based on smoothed estimators, such as those ones built with kernels. The effective calibration of 
the proposed test in this setting demands a bootstrap strategy. In addition, the issue of bandwidth 
selection, a common drawback of the tests based on kernel smoothing, has to be faced. This sub- 
section is devoted to provide a joint solution to these two problems. 



In order to design a resampling procedure for test calibration, samples from the directional-linear 
variable under the null hypothesis of independence must be drawn. Given that the null hypothesis 
is nonparametric, a smooth bootstrap procedure will be considered. That is, samples under Hq 
can be obtained from the product of the marginal kernel density estimators, fx. t h p fz;g p , where the 
notation h p and g p refers to pilot bandwidths. The simulation is easily performed with the following 
algorithm. 

Algorithm 1. Let {(Xj,Zj)}™ =1 be a random sample from a directional-linear variable (X, Z). 
I Sample independently two indexes i,j in 1, ... ,n. 
// Sample X* from a vM(X h l/hj). 
Ill Sample Z* from a J\f(Zj,g%). 
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Then, the pair (X.*,Z*) has density fx;h p fz;g p - 

It should be noticed that the resampling pilot bandwidths (h p ,g p ) cannot be used directly for es- 
timation purposes. The use of the same bandwidths for resampling and estimation leads to an 
inconsistent bootstrap procedure in the linear setting, where Cao (1993) suggested solving this 
problem by adjusting the degree of smoothness in resampling and estimation. Hence, two pairs 
of bandwidths are needed: (h p ,g p ) pilot bandwidths for resampling and (hb ,gb ) for estimation, 
where the subscript bo indicates that this pair is obtained by bootstrap bandwidth selection, as it 
will be explained later. 

The pilot bandwidth for the linear component, as suggested by Cao (1993), will be taken as: 

_ [ R{K") 1 ? 

where R(f") = J M f'z(%) 2 dx, H2{K) = f R K 2 (z) dz, being fz the marginal linear density. If the 
kernel function is a standard normal density, then R(K") = and ^(K) = 1. Since the true 
marginal density is not known, an approximate estimator for the integral of its third derivative must 
be obtained. This can be done by considering a reference normal model with standard deviation 
a, leading to R(f'z) = 16 ^| CT 7 ■ In practice, a robust estimator of a can be plugged-in the formula. 

It should be noted that g p = 0(n~?), in contrast with the optimal bandwidth minimizing the 

Asymptotic Mean Integrated Square Error (AMISE), which is 0(n~s) (Wand and Jones 1995). In 
that sense, the pilot bandwidth provides an extra smoothness with respect to the AMISE bandwidth. 



Unfortunately, there is no pilot bandwidth available for the directional case. To solve this problem, 
a possibility is to consider the AMISE bandwidth given in Garcfa-Portugues et al. (2012b) for the 
directional kernel density estimator (1) and settle a further degree of smoothness for the pilot band- 

width h p . The AMISE bandwidth for (1), which is 0{n 4+9 ), involves some quantities depending 
on the kernel and the underlying directional density, but with an explicit expression when consider- 
ing a von Mises kernel and the von Mises model as a reference density (see Garcfa-Portugues et al. 

(2012b) for details). Specifically, a plausible conjecture is that h p = Oyn 6 +<?) and therefore a 

_J 1_ 

possibility to select the pilot directional bandwidth is taking h p = /iamise • n 4 +<? 6 +<? , where /iamise 
is the AMISE bandwidth. 



Once the pilot bandwidths (h p , g p ) are obtained, appropriate estimation bandwidths for the directional- 
linear density must be also chosen. Considering the von Mises-normal kernel, the bootstrap ex- 
pression for the Mean Integrated Squared Error (MISE) of the directional-linear kernel density 
estimator (3) was derived by Garcfa-Portugues et al. (2012b): 

MISE* Vft) (h,g) = (c q (l/h 2 ) 2 C q (2// l 2 )- 1 2vr^n)~ 1 

+ n- 2 i' n [(i - n- l )** 2 (h) o n* 2 (g) - 2*i(/o o n*(g) + o n*] i n , 

where matrices ^f^(h) and fi*(g), a = 0, 1, 2 are given by 

_/ C q {l/hj) 2 \ _ ( f C q (l/h 2 )C q (l/h 2 p ) 2 e x ' x ^ h ' \ 

*° " ^,(||X l + X J ||/^J.. ' * W - [L q C q (\\x/h 2 + X t /h 2 \\) ^ (dx) ^/ 

_( f C q (l/h 2 ) 2 C q (l/h 2 p ) 2 \ 

n{k) ~ { k C q (I |x/tf + Xi/hjl |) C q (\\x/h 2 + X s /hj\\) ^ (dX/ ' 
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with 4>a a)g the normal density with standard deviation a a:9 = (ag 2 + 2g p r j 2 , a = 1, 2. The estimation 
bandwidths are obtained as the minimizers of the bootstrap MISE: 

(h bo , g bo ) = arg min MISE)^ (h, g) . 

To summarize, the pilot bandwidths (h p ,g p ) will be used for obtaining the estimation bandwidths 
(hb , gt,o), by minimizing the bootstrap MISE which has an exact expression (not requiring resam- 
pling). In addition, these pilot bandwidths will be also used in the following algorithm, in order 
to calibrate the test statistic distribution under the null hypothesis of independence between the 
directional and the linear components of the random variable. 

Algorithm 2. Let {(Xj, Zj)}™ =1 ^ e a random sample from a directional-linear variable (X, Z). 
I Compute the marginal pilot bandwidths h p and g p . 
II Obtain (h bo ,g bo ) = argmin M>0 MISE^ (h,g). 

III Compute the observed value ofT n from (5), with kernel density estimators taking bandwidths 
(hbo,9bo)- 

IV Bootstrap calibration. For b = 1, . . . , B: 

• Draw samples {(X* 6 , Z* 6 )}" =1 from fx-,h p fz;g p - 

• Compute T* b from (5), with kernel density estimators taking bandwidths (hb ,9bo)- 

V Approximate the p-value by ff{T n < T* b }/B, where if denotes the cardinal of the set. 

Both in Steps 3 and 4, instead of using the simple expression for T n given by (5), one might 
directly compute the original test statistic in its integral form. The same strategy could be used for 
computing the bootstrap bandwidths, where the direct definition of the bootstrap MISE could be 
used, or replaced by the previous formula in terms of matrices, available for the von Mises-normal 
kernel. In Step 4, sampling from fx;h p fz-,g p (under the Hq) is achieved using Algorithm 1. Note 
that the p-value is approximated by the proportion of bootstrap replicates of the test statistic that 
exceed its observed value in the sample. 

3.3 Some simulation results 

The finite sample performance of the test statistic and the proposed calibration algorithm is il- 
lustrated through a simulation study. Setting the null hypothesis of independence between the 
directional and the linear components, a family of distributions with a parameter controlling the 
departure from the independence has been considered. Specifically, the family of directional-linear 
densities in this study is given by: 

/(x,z)(x,*0 = (1 - S)f^(x)f z (z) + 5f A f\ vM (x,z), 

where the directional component /x will be taken as a von Mises model vM((O q , 1), 1) and the 
linear one is a standard normal. The null hypothesis of independence is satisfied for 5 = 0, and 
departures from this model are obtained by setting different values of 5. Finally, these departures 
are built by including f^ vM (x, z), a directional-linear density with non-independent marginals: 

, , , (z- (a(l + x / /i r )/2 + 6)\ 
/a^m(x, z) =<p I 1 /„ M (x; H, k). 

The previous density has a von Mises marginal in the directional component with mean fj, and 
concentration k, and the linear part is a normal density with standard deviation a but with mean 
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value conditioned on the directional part x. Intuitively, for a given x in Q q , the mean of the linear 
marginal depends on the position of x relative to the reference mean /j, r . If x = fi r , then the 
conditional mean of the linear component is (a + b), whereas if x = — fi r , the mean is b. See Figure 
3 for some examples in the circular-linear and the spherical-linear cases. 




Figure 3: Densities fm V M f° r the circular-linear and spherical-linear cases. Directional-linear data samples 
are shown as points in M. 2 (circular-linear case) and K 3 (spherical-linear case) as follows: distance from the 
origin represents the linear component, whereas direction of the point gives the directional component. For 
both cases, the von Mises marginal has mean /j, — (0 q ,l) (blue vectors) and concentration n = 1. The 
conditional mean of the linear component is 2 + x'/x r (corresponding to a — 2, b = 1), where fi r = (—1, q ) 
(red vectors) and its variance is <r 2 — 1. The blue points represent the directions sampled from the directional 
marginal and the red ones the final directional-linear sample. The red curve/surface is the conditional mean 
of Z given X = x. 

The algorithm presented in the previous section will be used in different simulated scenarios, under 
independence (namely, under Hq) and for different values of 5 (denoting these settings by Hs). 
Empirical size and power of the test will be explored, running 1000 Monte Carlo replicates and 
B = 1000 Bootstrap replicates. Density fj^\ v M will be used with parameters fi = (0 q , 1), k = 1, 
fi r = (—1, q ), a = 2, b = 1, and a 2 = 1. 

Percentages of rejections of the null hypothesis under Hq and Hs, with 5 = 0.10,0.20,0.30, are 
reported in Table 1, for the circular-linear and spherical-linear cases, with different significance 
levels and sample sizes. For the circular-linear case, results in Table 1 show that the empirical size 
is close to the nominal level for all the sample sizes considered. In addition, a satisfactory behavior 
in terms of power can be also concluded from the percentages of rejections under Hs- As expected, 
the number of rejections grows as long as the simulated model departs from the null hypothesis 
of independence, improving the results with increasing sample size. The results for the spherical- 
linear case are quite similar to the previous ones for the empirical size, but with lower power in 
comparison with the circular-linear scenario, something expected as a consequence of the curse of 
dimensionality. 
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Circular-linear Spherical-linear 



Models 


n 


a = 0.10 


a = 0.05 


a = 0.01 


a = 0.10 


a = 0.05 


a = 0.01 


H n nn 

-L-L U.UU 


200 


0.089 


0.053 


0.008 


0.089 


0.053 


0.008 




500 


0.091 


0.050 


0.015 


0.091 


0.050 


0.015 




1000 


0.095 


0.046 


0.012 


0.095 


0.046 


0.012 


/j (1 1(1 


200 


0.250 


0.163 


0.050 


0.194 


0.100 


0.032 




500 


0.634 


0.488 


0.203 


0.399 


0.279 


0.109 




1000 


0.934 


0.868 


0.650 


0.784 


0.666 


0.386 


/ / n on 


200 


0.767 


0.665 


0.392 


0.590 


0.464 


0.214 




500 


1.000 


0.996 


0.961 


0.971 


0.941 


0.815 




1000 


1.000 


1.000 


1.000 


1.000 


1.000 


1.000 


#0.30 


200 


0.987 


0.970 


0.878 


0.898 


0.819 


0.592 




500 


1.000 


1.000 


1.000 


1.000 


1.000 


0.998 




1000 


1.000 


1.000 


1.000 


1.000 


1.000 


1.000 



Table 1: Percentage of rejections for the test, for different significance levels a and sample sizes n. Circular- 
linear and spherical-linear cases. First row: null hypothesis of independence. Second to fourth row: alterna- 
tive models. 

4 Real data analysis 
4.1 Data description 

The original Portuguese fire atlas, covering the period from 1975 to 2005, is the longest annual 
and country wide cartographic fire database in Europe (Pereira and Santos 2003). Annual wild- 
fire maps were derived from Landsat data, which represents the world's longest and continuously 
acquired collection of moderate resolution land remote sensing data, providing a unique resource 
for those who work in forestry, mapping and global change research. For each year in the dataset, 
Landsat imagery covering Portugal's mainland was acquired for two time steps, pre- (early May) 
and post- (late September) fire season, thus providing a snapshot of the fires that occurred during 
the season. Annual fire perimeters were derived through a semi-automatic procedure that starts 
with supervised image classification, followed by manual editing of classification results. Minimum 
Mapping Unit (MMU), i.e., the size of the smallest fire mapped, changed according to available 
data. Between 1975 and 1983 (the MultiSpectral Scanner era), spatial resolution of satellite images 
is 80 meters and MMU of 35 hectares. From 1984 forward with data availability at spatial reso- 
lution of 30 meters (Thematic Mapper and Enhanced Thematic Mapper era) MMU is 5 hectares, 
allowing to map a larger number of smaller fires than in the 1975-1983 era. Below an MMU of 
approximately 5 hectares the burnt area classification errors increase substantially, and given the 
very skewed nature of fire size distribution, the 5 hectares threshold ensures that over 90% of total 
area actually burned is mapped. For consistency, and due to discrepancies in minimum mapping 
unit between 1975-1983 and 1985-2005, in this study only fire perimeters mapped in the latter were 
considered, which results in 26870 fire perimeters. 

This application is based on the watershed delineation proposed by Barros et al. (2012). In their 
work, watersheds were derived from the Space Shuttle digital terrain model (Farr et al. 2007) using 
the ArcGIS hydrologic toolbox (ESRI 2009). Minimum watershed size was interactively increased 
so that each watershed contained a minimum of 25 fire observations (see the cited work for more 
details). Fire perimeters were attributed to the watershed that contains its centroid. 

The orientation of fire perimeters and watersheds is determined by principal component analysis, 
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following the approach proposed by Luo (1998). Specifically, principal component analysis is applied 
to points that constitute the boundary of the object (fire or watershed), resulting an orientation 
given by the first principal component (PCI). The points of the boundary can be represented ei- 
ther in a bidimensional space defined by each vertex's latitude and longitude coordinates or in a 
tridimensional space, taking into account also the altitude. Then, the PCI corresponds to an axis 
that passes through the center of mass of the object and maximizes the variance of the projected 
vertices, represented in M 2 or in M 3 . The fact of computing the PCI also in M 3 aims for taking 
into account the variability of fires according to their slope, which, as the center plot of Figure 1 
shows, present marked differences between regions. Then, the orientation of the object is taken as 
the direction given by its PCI. 

It is important to notice that an orientation is not exactly a directional observation and that some 
conversion is needed for applying the directional-linear independence test. In the two-dimensional 
case, the orientation is an axial object (the orientation N/S is also S/N). These orientations can be 
encoded by an angular variable G [0, 7r), with period n, so 20 is an usual circular variable. Then, 
with this codification, the angles 0, vr/2, vr, 3vr/2 represent the E/W, NE/SW, N/S and NW/SE 
orientations, respectively. In the three-dimensional space, the orientation is coded by a pair of 
angles (©,<£) using spherical coordinates, where 9 G [0, 7r) plays the same role as the previous 
setting and <I> G [0, 7r/2) measures the inclination (only positive angles are considered as the slope 
of a certain angle u> equals the slope of —uj). Therefore, points with spherical coordinates (20, <!>), 
which lie on the upper semisphere, can be considered as a realization of a spherical variable. 

4.2 Results 

The main results of the independence test are collected in Table 4.2. The null hypothesis of inde- 
pendence between the wildfire orientation and the burnt area (in log scale, from now on) is rejected, 
either using orientations in M 2 or in M 3 (first column). The test is carried out using all the 26870 
observations for years 1985-2005 without regarding about the watersheds. The second column of 
Table 4.2 contains the p-values for the hypothesis of independence between the orientation of a 
watershed and the total burnt area of fires within the region. Independence is also rejected for a 
significance level of a = 0.05. 





Fires orientation 


Watersheds orientation 


PCI in R 2 


0.0000 


0.0420 


PCI in R 3 


0.0000 


0.0310 



Table 2: p-values for the independence test between the fires orientation and the burnt area (in log scale) 
and the watershed orientation and the total burnt area (in log scale) of the watershed. 

After concluding the presence of dependence between wildfire orientation and size, it is possible to 
carry out a spatial watershed analysis by applying the test to each of the watersheds, in order to 
detect if the presence of dependence is homogeneous or if it is only related with some particular 
areas. Figure 4 represents the maps of the p-values of the test applied to the observations of each 
watershed, using PCI in M 2 and in M 3 (left plot and right plot of Figure 4, respectively). The maps 
reveal the presence of 14 and 18 watersheds where the null hypothesis of independence is rejected 
with significance level a = 0.05, for the circular-linear and the spherical-linear cases, respectively. 
This shows that the presence of dependence between fire orientation and size is not homogeneous 
and it is localized at very particular watersheds. It is also interesting to note that the inclusion of 
the altitude coordinate in the calculus of the PCI leads to a richer detection of dependence between 
the wildfire orientation and size at the watershed level. 
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p-values 

■ [0.00,0.01) 

□ [0.01,0.05) 

□ [0.05,0.10) 

■ [0.10,1.00] 




Figure 4: p-values of the independence test for the first principal component PCI of the fire perimeter 
and the burnt area (in log scale), by watersheds. The left map represents the p-values for the circular- 
linear independence test with PCI in K 2 , i.e., computed from latitude and longitude coordinates of the fire 
perimeter. The right one represents the p-values for the spherical-linear independence test, with PCI in M 3 , 
computed from latitude, longitude and altitude coordinates of the fire perimeter. 

5 Discussion 

A nonparametric test for assessing independence between a directional and a linear component has 
been proposed, and its finite sample performance has been illustrated through a simulation study. 
Simulation results support a quite satisfactory behavior of the calibration algorithm for the test, 
although further study is required on bandwidth selection, both for resampling and estimation. 
Indeed, the resampling procedure could be improved by computing the pilot and bootstrap band- 
widths in each new sample, but at the expense of increasing the computational cost. 

The immediate application of the test to the total wildfire orientation and size data makes possi- 
ble the detection of dependence between these two variables, both for two-dimensional or three- 
dimensional orientation. The same conclusion holds for watershed orientation and total area burnt. 
A detailed study of each watershed allows for a more specific insight into the problem. In water- 
sheds where there is evidence of independence between fire size and fire orientation support that 
the event-based analysis (Barros et al. 2012) should yield results similar to those that would be 
expected of an area-based analysis. On the other hand, detection of dependence between fire size 
and orientation in watersheds with uniform orientation (Barros et al. 2012) highlights cases where 
there may be a mixture of orientations. In such cases an analysis taking fire size into account might 
find evidence of preferred directionality in fire perimeters. In watersheds where fire events show 
preferential orientation (non-uniform distribution) and there is dependence between size and orien- 
tation, fire orientation distributions are structured in relation to fire size, especially considering the 
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typically asymmetric nature of fire size distributions, dominated by a small number of very large 
events (Strauss et al. 1989). In these cases, an area-weighted analysis of fire perimeter orientation 
might lead different results than those found by Barros et al. (2012). 

From a technical perspective, it may be argued that the data may not be independent and identically 
distributed over space and time. Unfortunately, given the data gathering procedure (detailed at the 
beginning of Section 4) dependence patterns cannot be clearly identified. Accounting for temporal 
or spatial dependence directly in the directional-linear kernel estimator and in the testing procedure 
is still an open problem. 
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A Proof of Lemma 1 

Proof. The closed expression (just involving matrix computations) for T n is obtained by splitting 
the calculus into three addends: 

T n = I (/(X,Z);M(X, Z) - /jC;ft(x) fz-g{ Z )) w <?( dx ) dz 

= 1 (^X> (^) * (^) - JW«)/*m)W 



Lq A JK. \ -J j =1 

n n 



sst» w> m ■< (s*) * m ■< (^) 



= 1 " "9 

+ 



/ fX;h(x) 2 h;g(z) 2 ^g(rfx) dz 

JClqXB, 

(6) - (7) + (8). 



The first addend is: 

n n 



< 6 > ^SS U (rP) K (S 5 ) ' (rP) K m ^ 



14 



c q {i/h?y 



EE 



2y 



(Zi - Zj) 



1=1 J 

For the second addend, 
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Finally, the third addend is obtained as 

(8)= / f^ffz^zfw^dz 

= / fc-hW 2 Uq(dX-) / f Z ;g(z) 2 dz 

Jn q Jr 

n n 

EE 
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From the previous results and after applying some matrix algebra, it turns out that 



where: 
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